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Abstract 


A finite-difference computer program (COSIM) has been written which models the one- 
dimensional, diffusions transport associated with high-temperature oxidation and interdiffusion 
of overlay-coated substrates. The program predicts concentration profiles for up to three 
elements in the coating and substrate after various oxidation exposures. Surface recession due to 
solute loss is also predicted. Ternary cross terms and concentration-dependent diffusion 
coefficients are taken into account. The program also incorporates a previously-developed oxide 
growth and spalling model to simulate either isothermal or cyclic oxidation exposures. In 
addition to predicting concentration profiles after various oxidation exposures, the program can 
also be used to predict coating life based on a concentration dependent failure criterion (e.g., 
surface solute content drops to 2%). The computer code is written in FORTRAN and employs 
numerous subroutines to make the program flexible and easily modifiable to other coating 
oxidation problems. 


Introduction 

Many blades and vanes in gas turbine engines are coated with an aluminide or overlay 
coating to impart additional oxidation and corrosion protection to the component. These coatings 
provide protection by the selective oxidation of Al to form a compact and adherent Al 2 0 3 scale. 

Under isothermal conditions, the A1 2 0 3 scale thickens at a parabolic rate (scale thickness 
proportional to the square root of time). Diffusional transport within the coating supplies Al to 
the growing oxide scale at a rate consistent with the parabolic growth of the scale. Although Al 
is continually consumed from the coating, the rate is acceptably low such that the coating 
typically contains sufficient Al to easily provide protection during isothermal oxidation 
conditions. Although land-based turbines often operate in a nearly isothermal mode for long 
periods, aero gas turbines undergo thermal cycling with each flight. This thermal cycling, 
primarily to ambient temperatures when the engine is shut down, can cause cracking and partial 
spalling of the protective A1 2 0 3 scale. This spallation generally occurs randomly across flat 
component surfaces but is higher at edges and corners. The spalling may occur to the metal 
surface but more commonly occurs only in the outer layers of the oxide scale. However, loss of 
the oxide is not catastrophic since selective oxidation of Al continues when the component again 
reaches high temperatures such that the damaged scale will “heal” and continue to grow. 

Whenever scale spallation occurs due to thermal cycling during oxidation (i.e., cyclic oxidation), 
the oxide scale on the surface will, on average, be thinner than a scale grown isothermally for 
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equivalent hot exposure times. Since the rate of scale growth is inversely proportiona to the 
scale thickness, a consequence of this thinner oxide scale is that the average rate of Al 
consumption is greater during cyclic oxidation accompanied by scale spallation than during 
isothermal oxidation. Hence, Al is depleted from the coating at a higher rate during eye ic 
oxidation, the exact rate depending on the growth rate of the alumina scale and the amount ot 
oxide spalling. This paper will focus specifically on degradation of MCrAlY overlay coatings, 
where M stands for either Ni, Co or Fe. These coatings are often used on components m marine 
environments where hot corrosion protection is required, and in military applications where 
higher temperatures are often encountered. 

In addition to oxidative attack, overlay coatings are further degraded by interdiffusion with 
the substrate. Since the coating is, by nature, higher in Al content than the substrate, Al diffuses 
from the coating into the substrate and generally becomes unavailable to support the growth o 
the protective alumina scale. Likewise, the Cr concentration in the coating is also typically 
greater than that in the substrate resulting in the diffusion of Cr from the coating into the 
substrate. In contrast, Ni and other elements in the substrate diffuse into the coating. Althoug 
the Y in the alloy strongly affects adherence of the A1 2 0 3 scale during cyclic oxidation, it is at 
such low concentration in the coating (e.g., <0.2 at.%) that dilution by interdiffusion with the 
substrate is generally not considered. Schematic concentration profiles in the coating and 
substrate after a short exposure in an oxidizing environment are shown in Fig. 1 . When a coating 
is substantially depleted of Al during cyclic oxidation, it can no longer supply sufficient Al to re- 
heal the protective scale. Less protective oxides, such as NiO , FeO or CoO, can form signaling 
the end of the protective life of the coating. A critical Al content in the coating can be defined to 
indicate the useful life of the coating. This critical Al content could indicate the time at which 
less-protective oxides form on the surface, or an earlier time at which the coating could be 
stripped from the component and a new coating applied. Modeling the Al transport in the 
coating and substrate during cyclic oxidation allows the coating life to be predicted. The purpose 
of the present work was to develop a one-dimensional ternary diffusion model to predict the 
concentration profiles associated with the oxidation and interdiffusion of coated superalloys 
undergoing cyclic oxidation. This model, given the name COSIM for Coating Oxidation and 
Substrate Interdiffusion Model, employs finite-difference techniques embodied in a FORTRAN 
computer "program to provide numerical solutions to the appropriate diffusion equations. The 
computer code employs numerous subroutines to make the program flexible and easily 
modifiable to other high temperature coating oxidation problems. Although the computer 
program and discussion below are primarily in terms of a NiCrAl coating on a Ni-based 
substrate, other elements could be substituted for the Al and Cr. 

Ternary Diffusion Equations 

Because of the low Y content, MCrAlY overlay coatings can be approximated as ternary 
alloys with Al and Cr in a matrix of either Ni, Co or Fe. Ternary diffusion equations can be 
employed to simulate the Al and Cr transport associated with coating oxidation and 
coating/substrate interdiffusion. The choice of simulating Cr rather than Ni diffusion is 
inconsequential. In the ternary system, the concentration of the third component, being a 
dependent variable, is always determined by difference (i.e., C Ni =100 -C A ,-C Cr ). Since most 
superalloys are complex multicomponent, multiphase alloys, fully accounting for diffusional 
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interactions of the various superalloy components, as well as that of Y in the coating, on the A 
and Cr transport is beyond current capabilities. Hence, these potential interactions, other than 
those encountered in the ternary system, are ignored. In addition, NiCrAlY and CoCrAlY 
overlay coatings generally consist of two phases, a high-Al NiAl or CoAl P phase embedded in a 
Ni or Co solid solution y phase. This P phase is depleted as A1 is consumed as oxide and as A1 
diffuses into the substrate. It has previously been shown that this p phase is often completely 
dissolved well before the useful life of the coating has been reached (Ref. 1). It was also shown 
in this same reference that good agreement was achieved between measured and predicted 
concentration profiles in the coating at extended times after the p phase had been dissolved by 
assuming a single-phase coating for the entire oxidation exposure. Hence, in the current model 
development, the two-phase coating will be represented as a single phase. 


Fick’s laws describing diffusion in a ternary alloy are: 
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2nd Law Eq. 2 


where Jj is the flux of component j, C is the concentration of either A1 or Cr, D is one of the four 
ternary interdiffusion coefficients, or diffusivities, and X and t refer to distance and time, 
respectively. The first term on the right hand side (RHS) of Eq. 2 can be rewritten as: 
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Equation 3. a can be further expanded as: 


Eq. 3. a. 
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L j'k = Al,Cr Eq. 3.b. 


A similar equation exists for the second term on the RHS of Eq. 2, namely: 
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j,k = Al,Cr Eq. 3.c. 


Substituting Eqs. 3.b and 3.c into the two terms on the right hand side of Eq. 2 yields: 
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Eq. 4. 
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Initially, transport associated with oxidation and the interdiffusion of the coating and 
substrate may be viewed as independent problems. Hence, at short times, coating/substrate 
interdiffusion can be treated as interdiffusion between two semi-infinite materials at a location 
centered on the coating/substrate interface. Similarly, oxidation of the coating can be treated as 
oxidation of a semi-infinite material affecting only a thin region below the surface of the coating. 
These two regions, the “inner” diffusion zone resulting from coating/substrate interdiffusion and 
the “outer” diffusion zone resulting from transport associated with oxidation, are shown 
schematically in Fig 2. At longer times, the two diffusion zones will overlap and diffusion 
associated with the two regions must be considered together. A solution to Fick s second law (Eq 
2), whether analytical or numerical, requires initial and boundary conditions. The boundary 
condition at the oxide/coating interface is discussed in the following paragraphs whereas the 
initial conditions and other boundary conditions used by the COSIM model are discussed in a 
later section. 

For isothermal oxidation, the boundary condition at the oxide/coating interface, given as 
the rate of A1 consumption, is well defined and decreases uniformly with time (the time 
dependence of the rate is inversely proportionate to the square root of time). However, for cyclic 
oxidation, the rate of A1 consumption can vary in a non-uniform manner as the thickness of the 
oxide scale increases during the high temperature exposure but decreases on cooling as spallation 
occurs (Ref 2). In the present diffusion model, a separate model simulating oxide growth and 
spallation during cyclic oxidation has been adopted to provide the oxide/coating boundary 
condition. This oxide growth and spalling model (Ref 3), designated COSP by the authors, 
predicts the rate of A1 consumption during each cycle by continuously tracking the thickness of 
the oxide scale, accounting for growth during high-temperature exposures, and partial oxide loss 
on cooling. This rate of Al consumption ( J ox ) predicted by COSP is used as the boundary 
condition for the diffusion model. Hence, the supply of Al within the coating to the 
oxide/coating interface must equal J ox . However, since Al is consumed from the coating, the 
coating surface recedes due to the loss of matter. This recession (4) is given as (Ref 4). 

£ = - JW.J E<1 5 

where V A | is the partial molar volume of Al in the coating and J 0 x is the rate of Al consumption 
discussed above (which is also the flux of Al entering the oxide). J 0 % can also be considered as 
the Al flux to the left of the moving oxide/metal interface in Fig 2, or the Al flux away from the 
interface. Similarly, J A i is the flux in the metallic coating towards the interface. Because of the 
interface motion, J ox is greater than J At according to the relationship: 
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where £ and ^ refer to the left hand side and right hand side of the oxide/coating interface, 
respectively, as shown in Fig 2. The parameter a in Eq. 6 is given as: 

a ‘X'-v*'C Mfi ) Eq ' 7 

where C A , 0 is the A1 concentration in the coating at the oxide/coating interface (Fig 2). 
Obviously, for the hypothetical case for V At = 0 (the A1 atoms in the coating have no volume), 
a=\, the A1 flux to the interface, J A t, becomes equal to the A1 flux away from the interface, J ox , 
and the interface is stationary. The partial molar volume for A1 was assumed independent of 
concentration due to a lack of available data for most alloy systems of interest (e.g., the y, Ni 
solid solution phase inNiCrAl alloys). 


Hence, the boundary condition for A1 at the oxide/coating interface is given by Eq. 6 
whereby the rate of Al consumption due to oxide formation, predicted by the COSP oxide model, 
J ox , is equated to the supply, or flux, of Al to the interface within the coating, J At , while taking 
into account the motion of the interface through the parameter a. The rate of Al transport within 


the coating to the oxide/coating interface (i.e., the Al flux, J A i) is given by Eq. 1, stated as: 
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Although no Cr is assumed to enter the A1 2 0 3 scale, the surface recession requires the diffusion 
of Cr and Ni away from the oxide/coating interface into the coating. The boundary condition for 
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Eq. 9 


where J c , is the flux of Cr in the coating away from the interface and C f> 0 is the Cr concentration 
in the coating at this interface (Fig 2). In the computer program, the calculations for COSP have 
been contained in a subroutine to facilitate the incorporation of other oxide growth and spalling 
models (e.g., Ref 5). 


Finite-Difference (F-D) Method 

The initial step in the F-D technique is to establish a grid of equispaced nodes across the 
region of the material over which diffusion will occur (i.e., the diffusion zone). Each node has a 
specific concentration associated with it. Fick's laws (Eqs. 1 and 4) are replaced with F-D 
equivalents based on small differences in concentration, AC, distance, AX, and time, At. A 
solution yielding the concentration profile at some time t is derived by solving the appropriate F- 
D equivalents for small time increments (A t) in an iterative manner. These iterations, or time 
steps, are continued until the At increments sum to the desired time t. A portion of a diffusion 
zone at time t and the new concentrations at time t+At is shown schematically in Fig 3 (Ref 6). 
Because of the typically large number of repetitive and tedious calculations made each iteration, 
F-D solutions are ideally handled by a computer. 

The F-D equivalents to Fick's laws are based on Taylor series expansions (Refs 7,8). The 
F-D equivalent of Fick's 2nd law (Eq. 4) can be given by either an explicit or implicit 
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representation. The explicit form was used throughout this work. F-D expressions for both first 
and second order partial derivatives for concentration with respect to distance, given on the right 
hand side of Eq. 4, were given by first central difference equations, stated as: 

dC; _ C;>h- 1 ~ Cjjn - 1 £q JQ 

dX TAX 
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Eq. 11 


where j refers to either A1 or Cr, n refers to the node number and the superscript i refers to the 
current iteration at time t. These equations apply to all nodes n where nodes n-1 and n+1 exist. 
A first forward difference expression was used for dCjdt in the left hand side of Eq. 4, given 


as: 
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Eq. 12 


dt At 

where the superscript /+/ refers to the next iteration at time t+At. The time increment At for the 
explicit F-D method is limited by a stability criterion typically given as: 

At 
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Eq. 13. 


where D max is the appropriate diffusion coefficient which is discussed in a later section. 
Although At is initially very small, the grid expansion scheme described below allows the time 
increment each iteration to increase, generally allowing long term simulation of oxidation with a 
reasonable number of program iterations. 


The last terms in Eq. 4 to be discussed are the eight derivatives of the diffusion 
coefficients ( 8D/8C ). Obviously, the concentration dependence of the four ternary diffusivities 
is required to evaluate this expression. The concentration dependence of the diffusivities is input 
to the computer program in polynomial form. The computer program determines an abbreviated 
polynomial expression for each derivative using the polynomial coefficients. If some, or all of 
the ternary diffusivities are concentration independent, or if the concentration dependence is not 
known and cannot be input to the program, the value of the appropriate derivatives is zero. 
Hence, each of the terms on the right hand side of Eq. 4 have known values at time t and can be 
evaluated. Since At is determined from the stability criterion in Eq. 1 3 and values for C' „ are 
known, values for C£J from the left hand side of Eq. 4 for the new time (t + At ) can be 
calculated. 


For the boundary condition at the oxide/coating interface (Eqs. 8 and 9), no node exists in 
the oxide (i.e., no n-1 node exists, where n is located at the interface) so that the central 
difference formulas in Eqs. 10 and 1 1 cannot be used. Consequently, concentration gradients 
(dC i/dX ) at the interface were determined using second order forward difference equations. 
Hence, the F-D equivalent for the A1 concentration gradient in both Eqs. 8 and 9 at the interface 
was given as: 
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where the second subscript number refers to the node which is numbered sequentially from zero 
at the interface and the superscript refers to the current iteration at time t (i.e., C A/0 is the A1 

concentration in the coating at the oxide/coating interface, X = cf shown in Fig 2). A similar 
expression was used for the Cr gradient. The diffusivities for Eqs. 8 and 9 were evaluated for the 
A1 and Cr concentrations in the coating at the oxide/coating interface (i.e., C A i,o .Q>,o)- 
Substituting Eq. 14 into Eqs. 8 and 9 together with Eqs. 5-7 are sufficient to yield the interface 

concentrations C' AI0 ,C' Cr0 . 

COSIM Computer Program 

The COSIM program initially establishes separate diffusion zones for interdiffusion and 
oxidation. Typically, these diffusion zones are set to be very narrow (i.e., 0.1 micron) but are 
allowed to expand with increasing interdiffusion. When the diffusion zones eventually overlap 
within the coating, the two separate zones are combined and the simulation continues with a 
single zone across the coating and into the substrate. This approach allows the use of a 
reasonable number of nodes in a zone yet with a fine node spacing during the early times when 
the concentration gradients are steep. For greater accuracy, all variables used in the iterative 
calculations were defined as double precision. 

The starting width of the outer and inner zones, DXCOAT and DXSUB, and the number of 
nodes in each zone, NCOAT and NSUB, respectively, are parameters input to the program. The 
spacing between nodes, DELX1 for the outer zone and DELX2 for the inner zone (i.e., AX in the 
F-D equations), is equal to the width of the zone divided by the number of nodes minus one (i.e., 
NCOAT- 1 and NSUB-1, the minus one since both zones are bounded by nodes). The 
concentrations associated with each node for both diffusion zones are stored in arrays labeled 
‘AT and ‘Cr’. In the outer diffusion zone, the nodes are numbered sequentially starting with zero 
at the oxide/coating interface through NCOATH (equal to NCOAT-I). The NSUB nodes in the 
inner diffusion zone are numbered from NSUBL (low) to NSUBH (high) with NSUBL starting at 
a value of NCOAT+2. Central difference equations, such as Eqs. 10 and 1 1, cannot operate on 
the endpoints of a zone since concentrations at nodes n+1 and n-1 are required. Although 
forward and backward difference expressions can be used at these locations, a common technique 
is to add ancillary nodes to the zone to allow the continued use of central difference equations. 
Hence, an additional ancillary node was added to the inner end of the outer diffusion zone and 
an ancillary node was added to each end of the inner diffusion zone ( NSUBL -1 and NSUBH+ 1). 
The ancillary nodes maintain assigned constant concentration values. The node numbering 
scheme is schematically shown in Fig 4 for NCOAT=l (e.g., 7 nodes numbered 0 to 6) and 
NSUB=\ 1 (e.g., 1 1 nodes numbered 9 to 19). The ancillary node for the outer diffusion zone is 
shown as node 7 and the ancillary nodes for the inner diffusion zone are shown as nodes 8 and 
20. As shown, different values for DXCOAT or DXSUB , or for NCOAT or NSUB, may result in 
different values for DELX1 or DELX2. 
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Initial Conditions 

Initially, each of the nodes in the outer zone is assigned the concentration of the coating. 

The ancillary node (node 7 in Fig 4) maintains this coating concentration until the inner and 
outer diffusion zones overlap. These initial concentrations for the outer zone are shown in lg . 
An error function solution was used to assign the coating and substrate compositions to the nodes 
in order to provide a smooth transition between the coating and substrate compositions. The 
NSUBL to NSUBH nodes in the inner zone are assigned concentrations as: 

c, = c;™" + (Cj‘ h - C; M ' )*\{l + erf )) j = Al, Cr Eq. 1 5 

where C aH " , C"' h refer to the initial concentration of the coating and substrate, respectively and 
Xmod is a modified distance parameter ranging from -2 to +2 across the diffusion zone of width 
DXCOAT. This range in X mo d produces a smooth concentration gradient with end compositions 
within 0.5% of C‘""' ,C ■"* . The coating composition was assigned to the ancillary node at 
NSUBL-1 (node 8 in Fig 4) and the substrate composition was assigned to the ancillary node at 
NSUBH+1 (node 20 in Fig 4). Again, these ancillary nodes maintain these assigned 
concentrations until the diffusion zones overlap (i.e„ the ancillary nodes are not operated upon 
by Fick’s 2 nd law, as discussed below). The initial concentrations for the inner diffusion zone are 

also shown in Fig 4. 

Use of the error function solution to assign initial concentration values does not have a 
significant effect on later concentration values because of the small initial diffusion zone width 
and the small initial time increments. A value on the order of 0.1 microns was typically used for 
the initial diffusion zone width, DXSUB. Typically, 30 to 40 nodes are assigned to the zone so 
that the node spacing ( AX) is 0.0025 to 0.003 microns. Substituting this latter value into Eq. 1 3 
and using a typical value of 10' 10 cm 2 /s for D max yields a time increment At of only 0.000225 
seconds. Hence, hundreds to thousands of iterations will typically be performed before the first 
minute of simulated exposure time allowing ample iterations for the starting concentrations to 
adjust to satisfy Fick’s 2 nd law (Eqs. 2-4). Fortunately, when diffusion begins to affect the 
concentration at the ends of the zones, either at node NCOATH, NSUBL or NSUBH, the node 
spacings are allowed to expand (as shown in the following section) so that larger time increments 
can be used each iteration. Since the outer and inner zones may have different node spacings, the 
program determines the maximum diffusion coefficient (D max in Eq. 13) for either the coating or 
substrate composition and uses the smaller of the two node spacings to determine the smallest 
time increment, DELT (i.e., At in the F-D equations), according to Eq. 13. 

Surface Recession, Flexible Zones and Semi-Infinite Boundary Conditions 

The COSIM model utilizes a flexible grid technique to account for surface recession and 
to simulate the semi-infinite boundary conditions. As the outer surface recedes due to Al loss, 
the entire outer zone is shifted and the concentrations at each node are adjusted.^ The technique 
used to accomplish this shift and adjustment is referred to as a “Murray-Landis” (M-L) 
transformation (Ref 9). This transformation shifts the nodes an amount proportional to their 
position from the moving boundary (X=Q to maintain a uniform node spacing. The semi-infinite 
boundary condition for the outer diffusion zone can be stated as. 
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c. | = C, c '"" i = Al,Cr . 

This boundary condition is approximated by increasing the zone width whenever diffusion 
“significantly changes” the concentration at the node NCOATH. Each iteration, the 
concentration at this node is changed slightly in accordance with Fick s 2 law (Eq. 4). A 
“significant change” used in the COSIM program was taken as 0.005 at.%. Varying this value 
from 0.002 to 0.008 at.% changed the number of iterations required to reach a fixed time but had 
no significant effect on the predicted concentration profiles at short or long times (i.e., 1 and 
1000 hrs using the sample input data given in Appendix A). Hence, whenever the concentration 
of node NCOATH (either A1 or Cr) varied from the coating composition by 0.005 at.%, the zone 
width was expanded by DELX1 and all node positions and concentrations were adjusted 
according to the M-L transformation. To further illustrate this operation, A1 loss at the surface 
(node 0) and operation of Fick’s 2 nd law (Eq. 4) will eventually cause the initial A1 profile in the 
outer diffusion zone (Fig 4) to appear as shown schematically by the solid squares in Fig 5. 
Eventually, the A1 concentration at node 6 will decrease to a value 0.005 at.% below the 
concentration of the coating, C? . During this iteration, the total zone width, DXCOAT, will be 
increased by an amount equal to the node spacing such that. 

DXCOA T '= DXC OA T+DELX1 Eq. 16 

The number of nodes remains constant but the new node spacing is given as: 

DELXl'=DXCOATV(NCOAT-l) Eq. 17 

The node positions are shifted proportionately such that node 6 is shifted inward by DELX1 
while node 0 at the surface (X=g) undergoes no shift (Fig 5). The concentrations at each node 
are shifted in a similar manner with node 6’ taking the old position of the ancillary node 7 and 
being reassigned the concentration of the coating, Q" al (the concentration previously held by the 
ancillary node 7 at the same position). The ancillary node maintains the coating composition at 
the new position. 

The semi-infinite boundary conditions at either end of the inner diffusion zone are 
simulated in a like fashion. Hence, as diffusion changes the concentrations at nodes NSUBL or 
NSUBH in the inner diffusion zone by 0.005 at.%, the zone width is expanded by the node 
spacing, DELX2 and the node positions and concentrations are adjusted according to the M-L 
transformation. The expansion occurs at the end of the zone where the concentration change 
occurred such that the inner diffusion zone expands either into the coating or into the substrate 
for changes at node NSUBL or NSUBH, respectively. The ancillary nodes ( NSUBL-1 and 
NSUBH+ 1 ) maintain the substrate composition and are repositioned at the ends of the zones with 
the new node spacing. Whenever the zones are expanded or contracted, a new time increment 
DELT, is calculated according to Eq. 13. As discussed above, the smaller of the time increments 
calculated for the outer and inner diffusion zones is always used. Because of the concentration 
dependence of the diffusivities, one end of the inner diffusion zone may expand more than the 
other over the course of several thousand iterations. Although the inner and outer zones may 
appear very different with different widths, different numbers of nodes and different node 
spacings, both zones always operate with the same time increment such that the total exposure 
time for both the inner and outer diffusion zones is always identical. 
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The expansion of the inner and outer zones will eventually result in their impingement, or 
overlap, within the coating. From this time onward, the two diffusion problems, oxidation and 
coating/substrate interdiffusion, become coupled and can no longer be operated independently. 

At the time that the diffusion zones overlap, the COSIM model sums the current width of the 
inner and outer zones and redefines a single zone of equal width. Equidistant node spacings are 
also calculated using the combined number of nodes from the two zones (i.e., NCOAT+NSUB). 
The COSIM model then fits a natural cubic spline curve (Ref 10,1 1) through the concentration 
profiles from both zones and generates a single profile through the coating and into the substrate 
for both A1 and Cr at each of the new node positions. Figure 6a shows predicted A1 and Cr 
concentration profiles in the coating and substrate at the time (time=1.34 hrs) when the two 
diffusion zones overlap. The width of the inner diffusion zone is significantly larger with a 
larger number of nodes and a larger node spacing than that for the outer diffusion zone. Figure 
6b shows the same data after the spline interpolation. Note that the nodes are equally spaced 
across the entire diffusion zone. The spline subroutines are based on equations and code given in 
Ref 1 1 . Following the combination of the two diffusion zones, the COSIM model continues to 
simulate A1 transport to the surface and interdiffusion of the coating and substrate. Both 
boundary conditions, at the oxide/coating interface and in the substrate remain as before the 
combination. Loss of A1 continues to cause surface recession and a shrinking of the zone while 
diffusion in the substrate continues to result in zone expansion. Within the program, the 
transition from two zones to one is reflected in the value of the parameter ZONE. The value of 
ZONE changes from two to one after the diffusion zones are combined and the program operates 
on the single diffusion zone thereafter. Parameters are redefined so that the new, single diffusion 
zone utilizes parameters associated with the outer diffusion zone (e.g., DELX1, DELT1 , etc.) 

The model continues to simulate increasing oxidation exposure with each iteration of time. The 
program can print out concentration data at intermediate times (e.g., 50, 100, 500, 1000, 5000 
hrs) or at some predetermined failure condition (e.g., C AI 0 = 0 or 5 at .%) . Output is also 
written to files to ease plotting concentration/distance profiles. Figure 7 shows predicted 
concentration profiles after 1, 100 and 1000 hrs using the data given in Appendix A. Figure 8 
shows the time dependence of the surface concentrations AL(0), CR(0), rate of A1 consumption 
and total weight of A1 consumed which are also written to output files. Diffusivities for these 
examples were taken from reference 12. The value for the partial molar volume of A1 was taken 
from reference 13. 

Occasionally, oscillations in the values of the A1 concentration at the oxide/coating 
interface, AL(0), will occur the first few cycles. These oscillations can cause the computed value 
for the AL(0) to become negative or to exceed 100 at.%. It has also been observed that a very 
high initial flux from the COSP oxide growth and spalling model can also cause the computed 
value for AL(0) to be negative for hundreds of cycles. To protect against this physical 
impossibility, the program limits the values of AL(0) to be between 0 and 100 at.% (i.e., negative 
values are set equal to zero and values greater than 100 are set equal to 100). The consequence 
of these limits is that Eq. 6 is not satisfied during these iterations. However, monitoring the 
computed values for AL(0) typically show convergence to acceptable values at relatively short 
times (i.e., less than one hour when unreasonably high values for Kp were used with the example 
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input data given in Appendix A). Since these early oscillations could falsely trigger the test for 
AL(0) less than a critical A1 concentration (i.e., AL(0)<CRITAL) this test is not initiated in the 
program during the first two hours of oxidation exposure. No problems have been encountered 
using test cases using various values for CRITAL, including zero. 

Certain constraints on the time increment used each iteration ( DELT) are desirable when 
simulating certain oxidation conditions. For instance, a value for DELT much less than one hour 
might be desirable during cyclic oxidation with one hour cycles. In this case, a value of DELT of 
one or five minutes might be desired so that several iterations during oxide growth can be 
performed between periods of oxide spallation each cycle. In contrast, for isothermal oxidation 
for long periods, it might be preferable to have no constraint on the values for DELT in order to 
minimize the number of iterations to reach a solution. To provide this flexibility, a maximum 
value for DELT (MAXDT) may be input to the program. The default condition within the 
program is that no constraint be made on DELT beyond that given in Eq. 13. 

Program Application 

The program has been used to predict both concentration profiles and coating life after 
cyclic oxidation of overlay coated substrates (ref 1,14). The COSIM program has also been used 
to perform parametric studies to examine the effect of coating thickness and coating and 
substrate composition on coating life (ref 1). Recently, the program was modified to examine the 
coating life extension due to the presence of a perfect diffusion barrier between the coating and 
substrate (ref 14). 

The COSIM program was initially compiled and executed on mainframe computers and 
later revised for execution on a desktop PC. The current version of the program has been 
compiled with ANSI FORTRAN95. 

A description of the main program and each of the subroutines, as well as flowcharts, is given in 
reference 15. A description of each of the input parameters, an alphabetical list of all program 
variables, a sample input file and the corresponding output files are also given in this reference. 


References 

1 . J.A. Nesbitt and R. W. Heckel, Thin Solid Films , 119, 281,1 984. 

2. J.A. Nesbitt, “Diffusional Aspects of the High-Temperature Oxidation of Protective 
Coatings”, in Diffusion Analysis & Applications , Edited by A.D. Romig, Jr. And M.A. 
Dayananda, TMS, Warrendale, 1989, p. 307-324. 

3. C.E. Lowell, C.A. Barrett, R.W. Palmer, J.V. Auping, and H.B. Probst, Oxid. Met., 36, 81, 

1991. 

4. J.A. Nesbitt, J. Electrochem. Soc., 136 , 1518 (1989) 

5. K..W. Chan, Met. and Mat. Trans., 28 A, 411, 1997. 

6. J.A. Nesbitt, Oxid. Met., 44, 309, 1995. 

7. R.W. Hornbeck, Numerical Methods ( Quantum Publishers, New Y ork, 1975). 


February 16, 2001 


February 17, 2001 


11 



8. M.L. James, G.M. Smith, and J.C. Wolford, Analog and Digital Comp uter Methods, 
(International Textbook, Scranton, 1964) 

9. D. Murray and F. Landis, J. Heat Transfer, 81, 106, 1959. 

10. L.W. Johnson and R.D. Riess, Numerical Analysis , Addison-Wesley, Reading, MA (1977). 

1 1 . J.H. Mathews, Numerical Methods for Mathematics. Science and Engin eering, 2 nd Ed. 
Prentice Hall, Englewood Cliffs (1992). 

12. J.A. Nesbitt and R.W. Heckel, Met. Trans., 18A , 2075, 1987. 

13. J.A. Nesbitt, NASA TM 83738, Cleveland, OH, 1984. 

14. J.A. Nesbitt and Jih-Fen Lei, “Diffusion Barriers to Increase the Oxidative Life of Overlay 
Coatings” in Elevated Temperature Coatings: Science and Technology III, Edited by J.M. 
Hampikian and N.B. Dahotre, TMS, Warrendale, USA, 1999, p. 131-142. 

15. J.A. Nesbitt, NASA TM-2000-209271, Cleveland, OH, USA, 2000. 


Appendix A. 

Al concentration in the coating (at.%) 1 3.0 

Cr concentration in the coating (at.%) 20.0 

Al concentration in the substrate (at.%) 5.0 

Cr concentration in the substrate (at.%) 10.0 

Density of the coating (gm/cm 3 ) 7.754 

Partial molar volume of Al in the coating (cnT/mole) 7 . 1 

Coating thickness (microns) 100.0 

Number of nodes in the “outer” coating diffusion zone 30 
Number of nodes in the “inner” coating/substrate diffusion zone 40 
Initial width of the outer diffusion zone (microns) 0. 1 
Initial width of the inner diffusion zone (microns) 0.1 

Parabolic oxide growth parameter due to weight of oxygen in the oxide (mg /cm /hr) 
The hot cycle duration for each thermal cycle, in hours. 1 .0 
The spall parameter for the COSP spalling model (Ref. 3) 0.008 

The four ternary diffusion coefficients: 

D A uI(Cm C Cr )-[T229+(0. 0731 *C AI )+(-0. 0083 *C Cr )+(0. 01 01 *C A ?) + 

(0.000 16*Ccr 2 )]*10- 10 

D A ic r (C A i, C Cr ) =[0. 01 1 6+ (0. 0923 *C Al )+(-0 . 001 0 *C Cr )+(0. 0001 6 *C A , i 2 ) + 
(0.00001 7*C C r 2 )]* 10- 10 

D C rAi(C A i, C Cr ) —[ 0. 0766+(-0. 0153 *Cai) + (0. 0837*C Cr ) + (0. 00062 *C A ? ) + 
(-0.001 5*C Cr 2 )]* Hr 10 

DcrCr(C A i, Ccr) =[0. 783+(-0. 0123 *Cai) + ( 0. 0247*C Cr ) +(0. 00096 *C A i 2 ) + 
(-0.00057*C Cr 2 )]* lO' 10 
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Oxidation-Affected Interdiffusion-Affected 
Zone Zone 



Figure 1 . Schematic concentration profiles after oxidation exposure. Arrows indicate direction of atomic transport 
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Figure 2 Schematic concentration profiles after a short oxidation exposure. 
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Figure 3. Schematic diffusion zone for time t (squares) and time t+At (circles). The concentration subscript refers to the 
node number and the superscript refers to the time iteration. 
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Figure 5. Schematic Al concentration profile in the outer diffusion zone before and after the M-L transformation. 
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Figure 6a. Al and Cr concentration profiles at the time (Bme= 1.34 hrs) when the inner and outer diffusion zones overlap 
Input data for run given in Appendix A. (Alternate nodes hidden for clarity.) 





a 


Distance (microns) 

Figure 6b. Al and Cr concentration profiles at ttns-1 .34 hrs following the spline operation. Input data for run given in Appendix A 
(Alternate nodes hidden for clarity.) 
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Figure 7a. Al, Cr and Ni concentration profiles at time='].0 hr using the input data given in Appendix A. 
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Figure 7b. Al, Cr and Ni concentration profiles at f/'me=100.0 hrs using the input data given in Appendix A. 



U0pB^U33U03 Ifsl 


o to o in o m 

O 00 00 I s - I s - CO © 



in o m o in o 

<N (N — — 


(% m yz ) uo;^bu^u3duo3 .13 jo iv 


Figure 7c. Al, Cr and Ni concentration profiles at f/me=1 000.0 hrs using the input data given in Appendix A. 
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(File 0UT15) for the input data given in Appendix A. 





(jq/iiuD/Sui) sjB’a 


Time (hrs) 

Figure 8b. Time dependence of WMDOT (rate of Al consumption) from data written to unit 1 5 (File OUT 1 5) for the input 
data given in Appendix A. 
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